library(ncdf4)
## Flujo de radiación, onda larga
dlwrf <- nc_open("../../data/train/dlwrf_sfc_latlon_subset_19940101_20071231.nc")
## Flujo de radiación, onda corta
dswrf <- nc_open("../../data/train/dswrf_sfc_latlon_subset_19940101_20071231.nc")
stations <- read.csv("../../data/station_info.csv")
energy_train <- read.csv("../../data/train.csv")
Todos los predictores tienen el mismo formato:
3 variables:
Cada variable tiene 5 dimensiones:
print(attributes(dlwrf$var))
$names
[1] "intTime" "intValidTime" "Downward_Long-Wave_Rad_Flux"
print(attributes(dswrf$var))
$names
[1] "intTime" "intValidTime"
[3] "Downward_Short-Wave_Rad_Flux"
# Ver el formato
print(dlwrf)
File ../../data/train/dlwrf_sfc_latlon_subset_19940101_20071231.nc (NC_FORMAT_NETCDF4):
3 variables (excluding dimension variables):
int intTime[time] (Contiguous storage)
long_name: time as an integer (YYYYMMDDHH)
int intValidTime[fhour,time] (Contiguous storage)
long_name: valid time as an integer (YYYYDDMMHH)
float Downward_Long-Wave_Rad_Flux[lon,lat,fhour,ens,time] (Chunking: [6,3,1,4,1705]) (Compression: shuffle,level 4)
_FillValue: 9999
units: W m-2
long_name: Downward_Long-Wave_Rad_Flux_Average (Average for Mixed Intervals) @ surface
cell_methods: time: mean
GRIB_param_discipline: Meteorological_products
GRIB_param_category: Long-wave_Radiation
GRIB_param_name: Downward_long_wave_rad_flux
GRIB_generating_process_type: Forecast
GRIB_param_id: 2
GRIB_param_id: 0
GRIB_param_id: 5
GRIB_param_id: 192
GRIB_product_definition_template: 8
GRIB_product_definition_template_desc: Average, accumulation, extreme values or other statistically processed value at a horizontal level in a time interval
GRIB_level_type: 1
GRIB_level_type_name: surface
GRIB_interval_stat_type: Average
GRIB_VectorComponentFlag: easterlyNortherlyRelative
5 dimensions:
time Size:5113
long_name: Time
units: hours since 1800-01-01 00:00:00
axis: T
lat Size:9
long_name: Latitude
standard_name: latitude
units: degrees_north
actual_range: 31
actual_range: 39
axis: Y
lon Size:16
long_name: Longitude
standard_name: longitude
units: degrees_east
actual_range: 254
actual_range: 269
axis: X
ens Size:11
long_name: ensemble
standard_name: ensemble
axis: E
description: 0 is control, other values are perturbation numbers
fhour Size:5
long_name: forecast hour
units: hours
7 global attributes:
Conventions: CF-1.0
title: Subset of data from 2nd-generation multi-decadal ensemble reforecast generated from the NCEP Global Ensemble Forecast System, mimicking version operational at NCEP/EMC circa mid-2012.
institution: NOAA Earth System Research Laboratory (ESRL)
source: NCEP GFS v 9.01, T254L42. Control initial conditions from CFSRR. Perturbed initial conditions from ETR. Model error simulated with STTP.
references: http://www.esrl.noaa.gov/psd/forecasts/reforecast2/index.html
history: Subset created 2013-01-15 19:20:02 UTC
comment: Original dataset generated on DOE's supercomputers at Lawrence Berkeley Laboratory through ALCC/ASCR grant.
Cada variable tiene exactamente las mismas cinco dimensiones. Del resumen impreso arriba podemos rescatar los siguientes valores importantes:
| Dimensión | Tamaño | Valor Mínimo | Valor Máximo |
|---|---|---|---|
| Latitud | 9 | 31 | 39 |
| Longitud | 16 | 254 | 269 |
| Ensemble | 11 | 1 | 11 |
| fHour | 5 | 12 | 24 |
| Tiempo | 5113 | 1700586 | 1823256 |
Representan valores en la zona geográfica de interés.
En esta área, solo toman los siguientes valores:
Mapa de las estaciones meteorológicas y los sitios de medición
En estos datos, la longitud está en grados positivos desde el meridiano cero, llegando hasta 360°. Por otro lado, las coordenadas de las estaciones meteorológicas tienen este dato en valores que van de 0° a 180°, ya sea al Este u Oeste.
head(stations)
Esto será tomado en cuenta al etiquetar la longitud en los datos de entrenamiento.
# Primero obtener las variables de interés
dlwrf.data <- ncvar_get(dlwrf, "Downward_Long-Wave_Rad_Flux")
dswrf.data <- ncvar_get(dswrf, "Downward_Short-Wave_Rad_Flux")
# Etiquetar las dimensiones Latitud y Longitud
# Latitud: 9 [Índice 2]
# Longitud: 16 [Índice 1]
dim(dlwrf.data)
[1] 16 9 5 11 5113
# Valores obtenidos del resumen
lat.inicial = 31
lat.final = 39
# Se resta 360 para convertir de grados positivos a la representación que tienen
# las coordenadas de las estaciones.
lon.inicial = 254 - 360
lon.final = 269 - 360
dimnames(dlwrf.data)[[2]] <- seq(from=lat.inicial, to=lat.final)
dimnames(dlwrf.data)[[1]] <- seq(from=lon.inicial, to=lon.final)
dimnames(dswrf.data)[[2]] <- seq(from=lat.inicial, to=lat.final)
dimnames(dswrf.data)[[1]] <- seq(from=lon.inicial, to=lon.final)
# Muestra de como queda etiquetado
print(dlwrf.data[,,1,1,1])
31 32 33 34 35 36 37 38 39
-106 247.019 234.019 211.019 196.019 192.019 186.019 179.019 177.019 182.019
-105 241.019 233.019 217.019 203.019 197.019 187.019 181.019 190.019 190.019
-104 245.019 245.019 231.019 220.019 216.019 207.019 199.019 210.019 204.019
-103 256.019 247.019 232.019 229.019 227.019 220.019 214.019 218.019 213.019
-102 258.019 249.019 242.019 235.019 231.019 228.019 225.019 224.019 224.019
-101 269.019 255.019 251.019 243.019 238.019 234.019 234.019 232.019 232.019
-100 291.019 262.019 260.019 255.019 247.019 240.019 240.019 237.019 237.019
-99 327.019 280.019 266.019 261.019 251.019 245.019 244.019 240.019 240.019
-98 340.019 300.019 275.019 269.019 256.019 249.019 247.019 242.019 243.019
-97 349.019 330.019 290.019 277.019 268.019 259.019 249.019 243.019 246.019
-96 351.019 343.019 310.019 296.019 278.019 277.019 262.019 250.019 246.019
-95 337.019 337.019 321.019 314.019 292.019 279.019 270.019 261.019 255.019
-94 324.019 335.019 324.019 329.019 305.019 293.019 279.019 267.019 262.019
-93 336.019 336.019 331.019 330.019 314.019 310.019 302.019 282.019 266.019
-92 338.019 334.019 332.019 329.019 319.019 314.019 312.019 301.019 280.019
-91 323.019 328.019 326.019 322.019 319.019 316.019 311.019 310.019 297.019
De acuerdo con el glosario de la agencia que provee el conjunto de datos, un ensemble se define como:
Colección de modelos numéricos que muestran resultados ligeramente diferentes.
e1 <- dlwrf.data[,,,1,]
e2 <- dlwrf.data[,,,2,]
e3 <- dlwrf.data[,,,3,]
e4 <- dlwrf.data[,,,4,]
e5 <- dlwrf.data[,,,5,]
comparativo <- c(e1[1,1,1,1], e2[1,1,1,1], e3[1,1,1,1], e4[1,1,1,1], e5[1,1,1,1])
print(comparativo)
[1] 247.0190 246.0190 247.5191 248.0089 247.2623
El primer ensemble es el de control, y es el que será utilizado para este modelo.
dlwrf.data.e1 <- dlwrf.data[,,,1,]
dswrf.data.e1 <- dswrf.data[,,,1,]
dim(dlwrf.data.e1)
[1] 16 9 5 5113
Es la hora a la cual la estación meteorológica hace su pronóstico sobre la variable de interés. Esto sucede cinco veces en el día: Empieza a las 12, y se repite en intervalos de 3 horas, terminando a las 24.
dimnames(dlwrf.data.e1)[[3]] <- namedSeq("fhour_", 24, steps=3, start=12)
dimnames(dswrf.data.e1)[[3]] <- namedSeq("fhour_", 24, steps=3, start=12)
print(dlwrf.data.e1[1,1,,1])
fhour_12 fhour_15 fhour_18 fhour_21 fhour_24
247.019 239.000 242.000 262.019 263.019
Es la fecha en la que se tomó la medición de la variable. Esta columna corresponde directamente con una de las respuestas en los datos de entrenamiento.
Sus valores son obtenidos de la variable time contenida
en el conjunto de datos.
# Todas las variables manejan el mismo tiempo
time <- ncvar_get(dlwrf, "time")
# Diferente formato
head(time)
[1] 1700568 1700592 1700616 1700640 1700664 1700688
head(energy_train$Date)
[1] 19940101 19940102 19940103 19940104 19940105 19940106
Comparando la variable time con las fechas de los datos
de entrenamiento, se observa que su formato no coincide. Cada una está
representando su tiempo de forma distinta:
| Variable | Formato |
|---|---|
time |
Horas desde 1800-01-01 00:00:00 |
energy_train$Date |
Fecha formato YYYYMMDD |
Por lo tanto es necesario convertir nuestro tiempo en horas, a fecha.
dates <- as.Date(as.POSIXct(time*3600, origin='1800-01-01 00:00:00', 'UTC'))
dates <- gsub('-', '', dates)
# Además, energy_train$Date tiene sus fechas como ints
dates <- strtoi(dates)
head(dates)
[1] 19940101 19940102 19940103 19940104 19940105 19940106
head(energy_train$Date)
[1] 19940101 19940102 19940103 19940104 19940105 19940106
Con estos valores, se pueden etiquetar nuestros predictores:
dimnames(dlwrf.data.e1)[[4]] <- dates
dimnames(dswrf.data.e1)[[4]] <- dates
head(dlwrf.data.e1[1,1,1,])
19940101 19940102 19940103 19940104 19940105 19940106
247.0190 257.0190 231.8486 241.1632 242.1974 247.8309
En este primer experimento se hará un modelo solo para el sitio de generación ACME.
Para identificar los datos que son útiles para este sitio particular se utilizan las coordenadas:
print(names(stations))
[1] "stid" "nlat" "elon" "elev"
acme <- stations[stations["stid"] == "ACME"]
print(acme)
[1] "ACME" "34.80833" " -98.02325" " 397"
# Se redondea porque los archivos de los predictores tienen coordenadas enteras.
# Se tomará la medición más cercana.
acme.lat <- as.character(round(as.numeric(acme[2])))
acme.lon <- as.character(round(as.numeric(acme[3])))
print(c(acme.lat, acme.lon))
[1] "35" "-98"
Se determinó que las coordenadas más cercanas a ACME son: (35, -98).
Con estos datos extraemos las mediciones correspondientes.
acme.dlwrf <- dlwrf.data.e1[acme.lon, acme.lat,,]
acme.dswrf <- dswrf.data.e1[acme.lon, acme.lat,,]
# Para simplificar este primer experimento, se decidió promediar las predicciones
# a lo largo del día
acme.dlwrf.means <- colMeans(acme.dlwrf)
acme.dswrf.means <- colMeans(acme.dswrf)
acme.final <- data.frame(
Date = names(acme.dlwrf.means),
dlwrf = acme.dlwrf.means,
dswrf = acme.dswrf.means)
print(acme.final)
Ya que se tienen los predictores, solo hace falta conocer cuánta
energía se generó en ACME con esos valores de dlwrf y
dswrf.
energy.acme <- energy[c("Date", "ACME")]
colnames(energy.acme) <- c("Date", "ACME")
# Hay que tener el mismo formato para combinarlo con nuestros predictores
energy.acme <- transform(energy.acme, Date=as.character(Date))
print(energy.acme)
Y se combina con los predictores en un solo Dataframe para poder crear modelos.
energy.acme <- merge(energy.acme, acme.final, by="Date")
print(energy.acme)
Usando un solo predictor: dlwrf
Observaciones:
lm.fit <- lm(ACME ~ dlwrf, data=energy.acme)
plot(lm.fit)
summary(lm.fit)
Call:
lm(formula = ACME ~ dlwrf, data = energy.acme)
Residuals:
Min 1Q Median 3Q Max
-18039113 -4582784 656225 5697096 16272029
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1371158 591168 2.319 0.0204 *
dlwrf 47330 1777 26.639 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7375000 on 5111 degrees of freedom
Multiple R-squared: 0.1219, Adjusted R-squared: 0.1217
F-statistic: 709.6 on 1 and 5111 DF, p-value: < 2.2e-16
Observaciones
lm.fit2 <- lm(ACME ~ dlwrf + I(dlwrf^2), data=energy.acme)
plot(lm.fit2)
summary(lm.fit2)
Call:
lm(formula = ACME ~ dlwrf + I(dlwrf^2), data = energy.acme)
Residuals:
Min 1Q Median 3Q Max
-18877647 -4515620 436038 5530706 16824318
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.771e+07 3.452e+06 8.027 1.23e-15 ***
dlwrf -1.215e+05 2.188e+04 -5.554 2.93e-08 ***
I(dlwrf^2) 2.618e+02 3.381e+01 7.743 1.16e-14 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7333000 on 5110 degrees of freedom
Multiple R-squared: 0.1321, Adjusted R-squared: 0.1318
F-statistic: 388.9 on 2 and 5110 DF, p-value: < 2.2e-16
Verificamos que no haya correlación entre los datos
plot(energy.acme$dlwrf, energy.acme$dswrf)
lmm.fit <- lm(ACME ~ dlwrf + dswrf, data=energy.acme)
plot(lmm.fit)
summary(lmm.fit)
Call:
lm(formula = ACME ~ dlwrf + dswrf, data = energy.acme)
Residuals:
Min 1Q Median 3Q Max
-18874982 -1458088 1125228 2163414 21383542
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82952.2 301536.3 0.275 0.783
dlwrf -5292.1 1005.2 -5.265 1.46e-07 ***
dswrf 52945.4 438.8 120.663 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3759000 on 5110 degrees of freedom
Multiple R-squared: 0.7719, Adjusted R-squared: 0.7718
F-statistic: 8645 on 2 and 5110 DF, p-value: < 2.2e-16
lmm.fit2 <- lm(ACME ~ poly(dlwrf, 2) + poly(dswrf, 2), data=energy.acme)
plot(lmm.fit2)
summary(lmm.fit2)
Call:
lm(formula = ACME ~ poly(dlwrf, 2) + poly(dswrf, 2), data = energy.acme)
Residuals:
Min 1Q Median 3Q Max
-18628158 -1475913 955625 2254224 20281315
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16877462 51776 325.969 <2e-16 ***
poly(dlwrf, 2)1 -36602605 4386138 -8.345 <2e-16 ***
poly(dlwrf, 2)2 33081936 3738991 8.848 <2e-16 ***
poly(dswrf, 2)1 507944582 4166306 121.917 <2e-16 ***
poly(dswrf, 2)2 40333521 3982481 10.128 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3702000 on 5108 degrees of freedom
Multiple R-squared: 0.7788, Adjusted R-squared: 0.7787
F-statistic: 4497 on 4 and 5108 DF, p-value: < 2.2e-16